%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import warnings
warnings.filterwarnings('ignore')
plt.rcParams['figure.figsize'] = [10,7]

How to project 3D data to 2D maps

Projecting 3D data to 2D cannot be done without distortions

You can conserve either but never all:

• equal-angle/conformal: projected circles will stay circular
• equal-area: projected circles have same area
• equidistant: projected circles have same radii in conserving direction

ax = plt.axes(projection=ccrs.Orthographic())
ax.set_global()
ax.stock_img()
ax.gridlines()
ax.tissot()

<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb8e5453e20>

Cylindrical projections

• longitudes and latitudes are parallel
• longitudinal distances are equally-spaced and conserved

Angle conserving: Mercartor projection

ax = plt.axes(projection=ccrs.Mercator())
ax.stock_img()
ax.gridlines()
ax.tissot()

<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb8e6e57d90>

Area conserving: LambertCylindrical

ax = plt.axes(projection=ccrs.LambertCylindrical())
ax.stock_img()
ax.gridlines()
ax.tissot()

<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb8e6daefa0>

Distance conserving: PlateCarree projection

ax = plt.axes(projection=ccrs.PlateCarree())
ax.stock_img()
ax.gridlines()
ax.tissot()

<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb8e6d84340>

Good compromise without conservation: Miller projection

ax = plt.axes(projection=ccrs.Miller())
ax.stock_img()
ax.gridlines()
ax.tissot()

<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb8e6d512b0>

Pseudocylindrical/Elliptical projections

• typically are area conserving (but not all!)
• do not have parallel latitudes
• can have parallel longitudes

Area and distance conserving: Sinusoidal projection

ax = plt.axes(projection=ccrs.Sinusoidal())
ax.stock_img()
ax.gridlines()
ax.tissot()

<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb8e6d185e0>

Area conserving: Mollweide projection

• longitudes are parallel and unequally-spaced
• meridians are ellipses

ax = plt.axes(projection=ccrs.Mollweide())
ax.stock_img()
ax.gridlines()
ax.tissot()

<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb8e6c869a0>

Compromise with less distortion: Aitoff projection

• not equal-area, based on azimuthal equidistant projection
• longitudes are not parallel for less distortion at the edges
• edge is 2:1 ellipse

ax = plt.axes(projection=ccrs.Aitoff())
ax.stock_img()
ax.gridlines()
ax.tissot()

<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb8e6c51d00>

Area conserving with less distortion: Hammer projection

• equal-area, based on azimuthal equal-area projection
• longitudes are not parallel for less distortion at the edges
• edge is 2:1 ellipse

ax = plt.axes(projection=ccrs.Hammer())
ax.stock_img()
ax.gridlines()
ax.tissot()

<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb8e4427730>

Conclusion

Use area conserving elliptical projections like Mollweide or Hammer

# All you need to create projections is contained in the coordinate reference system
import cartopy.crs as ccrs

# In order to plot a map, one has to tell matplotlib which projection to use
ax = plt.axes(projection=ccrs.Mollweide())

# This is not an Axes object but a GeoAxes, a subclass of Axes representing a projection
print(type(ax))

# All normal matplotlib plotting methods work but have to include a transform in which format the data was defined
data = np.arange(4).reshape([2,2])
ax.imshow(data, transform=ccrs.PlateCarree())

<class 'cartopy.mpl.geoaxes.GeoAxes'>

<matplotlib.image.AxesImage at 0x7fb8e438e580>

# Be aware that cartopy uses the Earth ellipsoid as bases for projections

# You can redefine the globe to be a perfect sphere 
ax = plt.axes(projection=ccrs.Mollweide(globe=ccrs.Globe(ellipse='sphere')))

Arrays can be plotted normally with imshow

# Define projection
ax = plt.axes(projection=ccrs.Mollweide())

# Define grid where first index is latitudes and second index is longitudes
data = np.arange(72).reshape([12,6])

# As we defined a PlateCarree grid, we have to use the transform accordingly
# Imshow fits the grid assuming the data is regularly spaced from -180° to 180° longitude and -90° to 90° latitude
# However, imshow always flips the y-axis (origin='lower')
imap = ax.imshow(data, transform=ccrs.PlateCarree())

# Include colourbar
plt.colorbar(imap, location='bottom')

<matplotlib.colorbar.Colorbar at 0x7fb8e4306700>

Creating regular interpolated data of a spherical shell

# Load in grid data
# In this case this is 3d grid data on a regular cartesian grid
# However interpolation should work regardless of shape and size

density = np.load('rho.npz')['arr_0']
print(f'Density array has shape {density.shape}')
fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.set_xlabel('x [cells]')
ax1.set_ylabel('z [cells]')
ax1.imshow(np.log10(np.nansum(density, axis=1).T), origin='lower')
ax2.set_xlabel('x [cells]')
ax2.set_ylabel('y [cells]')
ax2.imshow(np.log10(np.nansum(density, axis=2)).T, origin='lower')
plt.tight_layout()

Density array has shape (159, 159, 159)

# Index grid
ind = np.indices(density.shape)

# Set origin to centre
centre = (np.array(density.shape)-1)/2.0
ind = (ind.T - centre).T

# Calculate radius
radius = np.sqrt(ind[0]**2 + ind[1]**2 + ind[2]**2)

# Cut out radial shell at given radius, here +-1 cell is enough -> flatten arrays
radcut = 75
mask = (radius <= radcut+1) & (radius >= radcut-1)
rad = radius[mask]
rho = density[mask]
x = ind[0][mask]
y = ind[1][mask]
z = ind[2][mask]

# Calculate longitude in degree
lon = np.rad2deg(np.arctan2(y, x))

# Calculate latitude (90 degree north) in degree
xyrad = np.sqrt(x**2 + y**2)
lat = np.rad2deg(np.arctan2(z, xyrad))

# Define equally-spaced interpolation points
step = 1.0
alon = np.arange(-180+step, 180, step)
alat = np.arange(-90+step, 90, step)
# Create grid from points
ilon, ilat, irad = np.meshgrid(alon, alat, radcut)

# Interpolate data
from scipy.interpolate import griddata
# Griddata interpolates the data with given points on the left to the new points on the right
irho = griddata((lon, lat, rad), rho, (ilon, ilat, irad))

# Define projection
ax = plt.axes(projection=ccrs.Mollweide())

# Plot interpolated data
# Centre is 0° longitude and 0° latitude
# Origin='lower' as imshow starts at top
# Cartopy is defined for geospatial data where the east is to the right
# In astrophysics we want to look from the 'inside-out' where east is to the west
# Therefore, we have to flip the longitude with extent
imap = ax.imshow(np.squeeze(np.log10(irho)), origin='lower', transform=ccrs.PlateCarree(), extent=(180,-180,-90,90))

# Include colourbar
plt.colorbar(imap, location='bottom', label=r'density [g cm$^{-3}$]')

<matplotlib.colorbar.Colorbar at 0x7fb8dd15d760>

# Interpolation does not know about cyclic data
ax = plt.axes(projection=ccrs.Mollweide()) 

# Plot shifted data, here centre is 180° longitude and 0° latitude
imap = ax.imshow(np.squeeze(np.log10(irho)), origin='lower', transform=ccrs.PlateCarree(), extent=(360,0,-90,90))
plt.colorbar(imap, location='bottom', label=r'density [g cm$^{-3}$]')

<matplotlib.colorbar.Colorbar at 0x7fb8dd06de20>

# Add cyclic data
large = np.where((lon>160) & (lon<180))
rho = np.append(rho, rho[large])
rad = np.append(rad, rad[large])
lat = np.append(lat, lat[large])
lon = np.append(lon, lon[large]-360)

small = np.where((lon<-160) & (lon>-180))
rho = np.append(rho, rho[small])
rad = np.append(rad, rad[small])
lat = np.append(lat, lat[small])
lon = np.append(lon, lon[small]+360)

# Redefine interpolation points as cyclic
step = 1.0
alon = np.arange(-180, 180, step)
alat = np.arange(-90+step, 90, step)
ilon, ilat, irad = np.meshgrid(alon, alat, radcut)
irho = griddata((lon, lat, rad), rho, (ilon, ilat, irad))

# Plot shifted data
ax = plt.axes(projection=ccrs.Mollweide()) 
imap = ax.imshow(np.squeeze(np.log10(irho)), origin='lower', transform=ccrs.PlateCarree(), extent=(360,0,-90,90))
plt.colorbar(imap, location='bottom', label=r'density [g cm$^{-3}$]')

<matplotlib.colorbar.Colorbar at 0x7fb8dcf5a250>

Plotting vector data

# Load velocity data
velazi = np.load('velazi.npz')['arr_0']
velpol = np.load('velpol.npz')['arr_0']

# Apply radial cut
vazi = velazi[mask]
vpol = velpol[mask]

# Add cyclic data
vazi = np.append(vazi, vazi[large])
vpol = np.append(vpol, vpol[large])
vazi = np.append(vazi, vazi[small])
vpol = np.append(vpol, vpol[small])

# Interpolate data
ivazi = griddata((lon, lat, rad), vazi, (ilon, ilat, irad))
ivpol = griddata((lon, lat, rad), vpol, (ilon, ilat, irad))

# Plot density map
ax = plt.axes(projection=ccrs.Mollweide()) 
imap = ax.imshow(np.squeeze(np.log10(irho)), origin='lower', transform=ccrs.PlateCarree(), extent=(180,-180,-90,90))
plt.colorbar(imap, location='bottom', label=r'density [g cm$^{-3}$]')

# Overplot quiver data
ax.quiver(-alon, alat, np.squeeze(ivazi), np.squeeze(ivpol), transform=ccrs.PlateCarree())

<matplotlib.quiver.Quiver at 0x7fb8dcee9a60>

Use regrid_shape to regrid quiver data

# Regridding is a feature of cartopy which interpolates vector data onto regular points in projection space
# The input can be either an integer which defines the minimum number of vectors in one direction or directly given by a tuple of ints

# Plot density map
ax = plt.axes(projection=ccrs.Mollweide()) 
imap = ax.imshow(np.squeeze(np.log10(irho)), origin='lower', transform=ccrs.PlateCarree(), extent=(180,-180,-90,90))
plt.colorbar(imap, location='bottom', label=r'density [g cm$^{-3}$]')

# Overplot quiver data
ax.quiver(-alon, alat, np.squeeze(ivazi), np.squeeze(ivpol), transform=ccrs.PlateCarree(), regrid_shape=44)

<matplotlib.quiver.Quiver at 0x7fb8dae7ce80>

Healpix

import healpy as hp

# Resolution is set by the parameter nside, which is generally a power of 2 
NSIDE = 2

# Get number of pixels
print(f'The number of pixels is: {hp.nside2npix(NSIDE)}')

# Healpy assumes the pixels are stored in a certain order in a flattened array
data = np.arange(hp.nside2npix(NSIDE))

# Ring order assumes horizontal rings starting at the North Pole
hp.mollview(data, title="Ring order")

The number of pixels is: 48

# Nested order assumes horizontal rings starting at the North Pole
hp.mollview(data, nest=True, title="Nested order")

# Healpix comes with very powerful functions

# Standard coordinates are co-latitude and longitude
# Co-latitide is 0 at the North Pole, pi/2 at the equator and pi at the South Pole
# Longitude starts from centre, increases to the west and goes from 0 to 2pi
print(f'Co-latitude and longitude of pixel 0 in radians: {hp.pix2ang(NSIDE, 0)}')

# Get longitude and latitude of pixels in degree, Healpy goes from 0° to 360°
print(f'Longitude and latitude of pixel 0 in degree: {hp.pix2ang(NSIDE, 0, lonlat=True)}')

# Get unit vector of pixels
print(f'Unit vector of pixel 0: {hp.pix2vec(NSIDE, 0)}')

# Get resolution of nside in radians or arcmin
print(f'Distance between pixels at NSIDE={NSIDE} is about {np.round(hp.nside2resol(NSIDE, arcmin=True)/60, 2)} degree.')

# Get size of pixels in square radians or square degree
print(f'Size of pixels at NSIDE={NSIDE} is about {np.round(hp.nside2pixarea(NSIDE, degrees=True), 2)} square degree.')

# Get pixel numbers of overlap with different shapes
print(f'Indices of pixels overlapping with disc with radius 0.1 in radians at z-vector: {hp.query_disc(NSIDE, [0,0,1], 0.1, inclusive=True)}')

Co-latitude and longitude of pixel 0 in radians: (0.4111378623223478, 0.7853981633974483)
Longitude and latitude of pixel 0 in degree: (45.0, 66.44353569089877)
Unit vector of pixel 0: (0.2825970826302196, 0.28259708263021954, 0.9166666666666666)
Distance between pixels at NSIDE=2 is about 29.32 degree.
Size of pixels at NSIDE=2 is about 859.44 square degree.
Indices of pixels overlapping with disc with radius 0.1 in radians at z-vector: [0 1 2 3]

Creating interpolated data of a spherical shell

# Increase the resolution to a similar one we used with cartopy
NSIDE = 64

# Define interpolation points, Healpy goes from 0° to 360°
ilon, ilat = hp.pix2ang(NSIDE, np.arange(hp.nside2npix(NSIDE)), lonlat=True)

# Shift data longitudes to be defined from 0° to 360°
slon = lon + 180

# Interpolate data
irho = griddata((slon, lat, rad), rho, (ilon, ilat, radcut))

# Plot data where we have to compensate for the 180° shift
hp.mollview(np.log10(irho), rot=(180,0,0), unit=r'density [g cm$^{-3}$]', title=None)

359.296875

Creating spherical integrations - how to distribute cell content

# Directly using cell centre gives bad results

# Index grid
ind = np.indices(density.shape)
# Set origin to centre
centre = (np.array(density.shape)-1)/2.0
ind = (ind.T - centre).T
# Calculate radius
radius = np.sqrt(ind[0]**2 + ind[1]**2 + ind[2]**2)

# Mask region
mask = (radius >= 1.0) & (radius<70.0)
rho = density[mask]
# Get pixel index of vectors
pix = hp.vec2pix(NSIDE, ind[0][mask], ind[1][mask], ind[2][mask])

# Sum up density in pixel
rhosum = np.zeros(hp.nside2npix(NSIDE))
rho = np.nan_to_num(rho)
np.add.at(rhosum, pix, rho)

# Plot result
hp.mollview(np.log10(rhosum), title='No distribution')

# Rough correction by spreading data to nearest neighbors, this does not work for all distances

# Sum up density in pixel
rhosum = np.zeros(hp.nside2npix(NSIDE))
rho = np.nan_to_num(rho)
np.add.at(rhosum, pix, rho/9)

# Add smoothed data by distributing cell content also onto next eight neighbors 
pixneigh = hp.pixelfunc.get_all_neighbours(NSIDE,pix)
for i in range(8):
    np.add.at(rhosum, pixneigh[i], rho/9)
    
# Plot result
hp.mollview(np.log10(rhosum), title='Simple smoothing onto eight neighbors')

# Better correction by getting the overlapping pixels with the cell

# Calculate angular extent of cells on the map in radians
# Here I estimate the mean angular extent of the regular cartesian grid to be sqrt(2) as we have all viewing angles
ext = np.arctan2(np.sqrt(2)/2, radius)

# Convert cell indices to vectors, flatten arrays
ext = ext[mask]
vec = np.array([ind[0][mask], ind[1][mask], ind[2][mask]]).T

# Get overlapping pixels and add data
rhosum = np.zeros(hp.nside2npix(NSIDE))
for i, ivec in enumerate(vec):
    ioverlap = hp.query_disc(NSIDE, ivec, ext[i], inclusive=True)
    np.add.at(rhosum, ioverlap, rho[i]/len(ioverlap))

# Plot result
hp.mollview(np.log10(rhosum), title='Simple pixel overlap')

# Further correction by weighting overlapping pixels with distance

# Calculate angular extent of cells on the map in radians
ext = np.arctan2(np.sqrt(2)/2, radius)

# Convert cell indices to vectors, flatten arrays
ext = ext[mask]
unitvec = (np.array([ind[0][mask], ind[1][mask], ind[2][mask]])/radius[mask]).T
ang = np.array(hp.vec2ang(vec, lonlat=True))

# Get overlapping pixels and add data
rhosum = np.zeros(hp.nside2npix(NSIDE))
for i, ivec in enumerate(unitvec):
    # Get indices of overlapping pixels
    ioverlap = hp.query_disc(NSIDE, ivec, ext[i], inclusive=True)
    # Calculate distance of overlapping pixel centre to vector in radians
    dist = np.arccos(np.dot(ivec, hp.pix2vec(NSIDE, ioverlap)))
    # Calculate weight as some function of the inverse of distance
    weight = (1/dist**0.4/np.sum(1/dist**0.4))
    np.add.at(rhosum, ioverlap, rho[i]*weight)

# Plot result
hp.mollview(np.log10(rhosum), title='Weighted pixel overlap')